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1 INTRODUCTION 



Young protostars can be highly variable (e.g., iHerbid 
Il983). The best known members of this club are the Fu 



Orionis objects, which are young stars showing very large 
(fact ors of hundreds) and abrupt increases in their bright- 
ness (|Hartmann fc Kenvon|[l996l ). Observations of these 
objects are best understood as variability in the accre- 
tion rate onto the growing protostar from a base accre- 
tion rate state of ~ 10~^ M0 yr~^ to a high of up to 
10~* Mq yr~^. Although detailed statistics is still lack- 
ing, the duration of the outbursts is believed to be in the 
range of tens to hundreds of years. Clear spectroscopic 
signatures implicate an accretion disc in the inner 1 
AU as the driver of this variable accretion rate and lu- 
mino sity output (|Zhu et al.l 1200 TI : lEisner fc Hillenbrandl 
I2OIII ). The rise times of the outbursts are very short, 
from a year to ~ ten years. FU Ori phenomenon is not 
just a fine detail of young stars evolution: the statistics of 
FU Ori objects suggests that an average proto-star may 
experience 10-20 of such outbursts, which may still be 
considered a lower limit on the import ance of the FU Ori 
outbursts (jHartmann fc Kenvon|[l99^ . It is not impossi- 



ABSTRACT 

It has been recently proposed that giant protoplanets migrating inward through 
the disc more rapidly than they contract could be tidally disrupted when they 
fill their Roche lobes '-^0.1 AU away from their parent protostars. Here we con- 
sider the process of mass and angular momentum exchange between the tidally 
disrupted planet and the surrounding disc in detail. We find that the planet's 
adiabatic mass-radius relation and its ability to open a deep gap in the disc 
determine whether the disruption proceeds as a sudden runaway or a balanced 
quasi-static process. In the latter case the planet feeds the inner disc through 
its Lagrangian LI point like a secondary star in a stellar binary system. As the 
planet loses mass it gains specific angular momentum and normally migrates 
in the outward direction until the gap closes. Numerical experiments show that 
planet disruption outbursts are preceded by long "quiescent" periods during 
which the disc inward of the planet is empty. The hole in the disc is created 
when the planet opens a deep gap, letting the inner disc to drain onto the star 
while keeping the outer one stalled behind the planet. We find that the mass- 
losing planet embedded in a realistic protoplanetary disc spawns an extremely 
rich set of variability patterns. In a subset of parameter space, there is a limit 
cycle behaviour caused by non-linear interaction between the planet mass loss 
and the disc hydrogen ionisation instability. We suggest that tidal disruptions 
of young massive planets near their stars may be responsible for the observed 
variability of young accreting protostars such as FU Ori, EXor and T Tauri stars 
in general. 

ble that proto-stars accrete most of its mass during these 
relatively short but very intense growth episodes. 

A rather natural model for FU Ori variables is a 
viscously and thermally unstable accretion disc that un- 
dergoes periodic switching be tween the stable quiescent 
and the outburst states (e.g., iBell fc Linlll994h . Hydro- 
gen is mainly neutral in the quiescent state and is mainly 
ionised in the outburst state, partially explaining the 
large difference in the accretion rates. In detail, how- 
ever, external triggering of the outbursts may still be 
desirable (e.g.. iBonnell fc Bastienlll992l : iBell et all [19951 : 
iLodato fc Clarke! |2004 ). as the theoretically produced 
outburst rise time scales are longer than ~ 10 yrs. 

The key defici ency of the pure thermal disc insta- 
bility model (e.g., iBeU fc LinI 1 1994 ) is in the fact that 
the unstable inner disc region is relati vely small, e.g., 
7? ~ 20i?Q ~ 0.1 AU (|Zhu et al.|[2007l ). This produces 
outbursts that are too short, unless one uses very small 
values for t he viscosity a-paramet er, e.g., as small as 
a = 10 ""^ (iLodato fc Clark j l2004l. This seems to be 
unlikely (|King et al.i.2007. : Zhuet al.ll2009l '). In addition, 
spectr al modeling of th e inner disc of the FU Ori vari- 
ables (|Zhu et al.r[2007l ), and their NIR interferometry 
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lEisner fc Hillenbrandl (|201ll ). require the hot disc region 
to be greater than R ~ 0.5 — 1 AU, which is clearly 
too large for th e pure thermal disc instability models 
(jZhu et al.ll2009l ). The small extent of the unstable region 
translates into too small a disc mass participating in the 
instability (since the mass is set by the crit ical surface 
densities at which the instability operates; iBeU fc LinI 
I1994I I. The disc instability model thus may be said to 
suffer from the "mass and time-scale deficit problem" as 
there is not enough mass in the disc to explain the ob- 
served long and bright FU Ori outbursts. 

A very promising way forward to enlarging the 
unstable region, and thus solving the deficit problem, 
was proposed by lArmitage et al] |2003). These authors 
pointed out that the intermediate disc region, R ~ 0.1 — 2 
AU is expected to have a layered structure (jGammid 
1 19961 ). The disc midplane region is shielded by the up- 
per layers from the ionisation sources, and is likely to 
be too co ld to support the Magn eto-Rotational Instabil- 
ity (e.g., iBalbus fc Hawlevlll998l ) that would otherwise 
yield an efflcient angular momentum transport. Accre- 
tion thus proceeds only inside the upper ionised "skin" of 
the disc. Bulk of the material is underneath the skin, too 
cold to flow radially, settling into what is called "the dead 
zone" . Over time, and in the right range of external mass 
supply rates, the dead zone becomes very massive, e.g., 
~ 0.1 M0, which eventually triggers a grav itational insta- 
bility . The associated gravito-turbulence (|Lin fc Pringld 
ll987l : lRice et al.ll2005l) initiate s a ra pid angular momen- 
tum transfer. Armitage et al.l (|200lh obtained relatively 



bright, M ~ 10~ M0 yr~ , outbursts lasting for as long 
as 10^ — 10^ yrs. The cu rrent consensus is that the model 
of lArmitage et al] (|200lh satisfies observational and phys- 
ical constr aints on FU Or i outbursts better than any 
other (e.g.. lZhu et~al]|2009l l. 

Here we propose another solution to the inner disc 
mass deficit problem, which relies on an entirely different 
mass reservoir. We suggest that the material accreted 
onto the star during the FU Ori outbursts comes from 
giant young planets that fill their Roche lobes inside 
the inner 1 AU from the star. The planets are assumed 
to be pushed inward closer to the star by gravitational 
torques of the protoplanetary disc. The significance of 
young planets as serious contenders on the role of mass 
donors for FU Ori outbursts is easily appreciated. The 
mass of gaseous proto-planets may be as high as tens of 
Jupiter masses (Mj), much higher than the likely disc 
mass in the inner fractions of an AU. If the material of 
such a planet is dumped into the inner disc and then ac- 
creted onto the star at an accretion rate of ~ 10~^ M0 
yr~^, then there is enough gas for a ~ O(IOO) yrs-long 
accretion outburst. 

Our work is based on ideas of iNavakshinI (|2011bh 
who pointed out that the early {t ^ 1 Myr) stages of the 
thermal relaxation of the planet are highly uncertain and 
model dependent (see a fuller discussion of this in ii3.4|l . 
The "low density start" models of such planets are suffi- 
ciently fluffy to be tidally disrupted inside a fraction of an 
AU. We sha ll also find important co nnections to the ear- 
lier work of iLodato fc Clarke! (|2004l ') who showed that a 
massive planet in the inner disc region can trigger rapid 
rise outbursts, as observed. Finally, our calculations of 



young protoplanet disruptions are complimentary to the 
"outburst acc retion mode" simulations of larg e (10—1000 
AU) discs by IVorobvov fc Basul (|2005l . l20od ). These au- 
thors have found that massive gaseous clumps born in 
their discs at ~ 0(100) AU distances migrate inward 
quickly, being accreted by their inner boundary condi- 
tion. 

The purpose of our paper is to investigate the pro- 
cess of tidal disruption of compact proto-planets in the 
inner ~ 0.1 AU of the star in unison with a more self- 
consistent modeling of the protostellar disc evolution to 
probe the proposed link of young protoplanets to the 
FU O ri outbursts. To model th e disc, we use the s tan- 
dard (|Shakura fc SunvaevI Il973l : iFrank et al] I2OO2I ') ID 
time-dependent viscous disc evolution equations in the 
presence of a massive embedded satellite. An approxi- 
mate but well understood and tested type II migration 
model is used to communicate t he gravitational torques 
between the planet and the disc (|Lin fc Papaloizoull 19861 : 
iLodato fc Clarke! |2004 ). In our modeling of the mass 
loss, we are aided by the well known results from semi- 
detached stellar binaries (e.g., [Rit tcr 1988). 

The structure of the paper and main results are the 
following. In ^ we build a simplified analytical theory 
for the mass loss and radial motion of the planet em- 
bedded in a proto-planetary disc. We show that due to 
conservation of angular momentum, a mass-transferring 
planet may migrate outward. Repeating the well known 
stellar binary results, we also find that the mass transfer 
rate and the migration pattern can be a quasi-steady or 
run-away process that may end in a complete (or nearly 
so) destruction of the planet. We introduce our numerical 
methods in ^ In |4] we present numerical experiments 
in which the outer disc is replaced by a constant external 
torque parameter. This allows us to test our simplified 
theory in a "clean" controlled setting. In !j5] simulations 
with a self-consistent disc evolution and torque on the 
planet are presented. A range of variability patterns is 
found there, and further expanded on in ^6] The effective 
temperature profiles of our disc models are discussed in 
^ In §8 we summarise the main results of our paper. 



2 ANALYTICAL CONSIDERATIONS 

The numerical scheme presented in |4]is well suited for 
studying tidal disruptions of planets in arbitrarily ini- 
tialised time-dependent gas discs. Here we attempt to 
build a simple analytical understanding of what to ex- 
pect from numerical simulations later on. We start with 
the simplest possible situations and gradually add differ- 
ent physical effects. In this section, the outer accretion 
disc is only considered as a source of a given gravita- 
tional torque on the planet (te, cf. §2.3) affecting its ra- 
dial motion. The inner disc interaction with the planet 
is described by the parameter fa introduced shortly be- 
low. We emphasise that re and fa are only defined for 
convenience of our analytical study (this section); these 
parameters are replaced by a self-consistent treatment of 
the disc-planet torques in the time-dependent numerical 
calculation. 
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2.1 Outward radial migration driven by 
planet's mass loss 

A planet embedded in an accretion disc is a sub- 
ject to gravitational torques fro m the surrounding discs 
that in general push it inwa rd (|Lin fc Papaloizou[|l979l : 
iGoldreich fc Tremaind Il98d ). Here we point out that 
a mass- losing planet can push itself outward. Opera- 
tionally, the effect is due to a different balance between 
the inner and the outer disc torques than that in the 
standard no-mass-loss case; however, the effect is best 
understood simply on the basis of angular momentum 
conservation in the star-planet system. 

Consider a planet of mass Mp orbiting the star of 
mass M, Mp at a distance a in a circular orbit. The 
planet's orbital angular momentum is given by 

Jp = MpQaa = Mp {GM^af^ . (1) 

If the planet fills its Roche lobe an d transfers mass in ward 
through the Lagrangian LI point (|Frank et al.ll2002l ') . the 
system is analogous to a semi-detached stellar binary sys- 
tem with an extreme mass ratio. In the case of conser- 
vative mass transfer, we have dM^/dt — —dMp/dt > 0. 
We neglect the spin angular momentum of the planet as 
it is small compared to the orbital angular momentum. 
We first assume that the star's spin does not change as a 
result of gas accretion onto it, in which case the planet's 
angular momentum is conserved during the mass transfer 
process. Therefore, 



(2) 



d In a d In M„ 
= -2 . 

dt dt 

This demonstrates that the planet migrates outward as 
it loses mass. This effect is well known in stellar binary 
evolution where a conservative mass exchange between 
the two components leads to widening of the orbit if the 
secondary is less massive than the primary, and shrinking 
of the orbit if mass is lost by the more massive star ijRitteil 
Il988l ). 

In a more general case, accretion of mass onto the 
star may spin up the star. Alternatively, the star could 
be slowing down its rotation by passing angular mo- 
mentum to the inner disc via magnetospheric torques 
(ICollier Cameron fc CampbeiilllQQa : [Armitage fc Clarke! 
Il996l ). Also, as we shall see later, the angular momen- 
tum of the gas accreting onto the star may flow past 
the planet's orbit if the gap, which is generally opened 
in the disc by a massive planet, is closed. If the viscous 
time of the inner disc is long compared with the planet's 
mass loss timescale, then we should also account in equa- 
tion [2] for the angular momentum stored in the disc. To 
continue with our analytical arguments below, we find it 
convenient to introduce a free parameter describing the 
efficiency of the outward migration due to mass transfer, 
fa > 0: 



din a d In Mp 

= —^Ja 



(3) 



dt ' dt 

If the star is a sink of angular momentum, fa < 1, but if 
it gives the angular momentum to the disc then fa may 
be larger than unity. 

We emphasise that, while introduced for convenience 
in the analytical considerations described here, fa is not a 



free parameter of our numerical simulations (see Section 
3 below), since the outward torque on the planet due 
to the inner disc is calculated self-consistently. We only 
use fa in our analytical model below, and also for an 
interpretation of the numerical results in cases where we 
can make a good guess on the appropriate value of fa. 



2.2 Planet's mass loss rate 

Mass transfer from a Roche lobe-filling secondary to the 
primary has been a subject of numerous papers in the 
context of cataclysmic binaries secular evolution (e.g., 
lRitteij[l988l ). The main result of these studies is that the 
mass loss rate is a very strong function of Ar = — r^f , 
where 



r-H 



: a(Mp/3M,)^/^ 



(4) 



is the Hills radius of the planet, which we set approx- 
imately equal to its Roche lobe radius, and Vp is the 
planet's radius. The mass loss rate is very small for neg- 
ative Ar, increasing rapidly for positive Ar. 

These basic facts are in fact sufficient for the analyt- 
ical treatment below, but we find it timely to introduce 
our numerical mass transfer scheme here as well, to em- 
phasise the importance of the difference Ar once more. 
Before the Roche lobe of the secondary is filled, Ar < 0, 
the mass transfer rate has an exponential form reflecting 
the exponential decrease of d ensity with height in the 
stellar atmosphere (cf. eq. 9 of lRitteJll988l l: 



dM, 



p _ PphCphr_f//ip 



dt 



»l/2 



exp 



Ar 



(5) 



where hp is the scale height of the planet's atmosphere, 
Pph and Cph are the photosphere's density and sound 
speed, respectively. 

When the Roche lobe of the planet is filled, and 
Ar > 0, the mass transfer rat e can be calculated i n 
a similar way (see section 1 of IPringle fc Wadd [l985l ') . 
and the result turns out to be dependent on whether 
the outer layers of the sta r are convective or radia- 
tive (|Rappaport et all Il983l : iKolb fc Ritteil Il992l ). For 
an isothermal atmosphere, dMp/dt oc exp[Ar /hp], with 
hp =const, whereas for a convective star with polytropic 
index 7, the result is 



dM, 



dt 



P ^ _^j.(37-l)/(2{7-l)) 



if Ar > 



(6) 



(equation (1.27) in lPringle fc Wadd[l985h . Importantly, 
for any value of 7 > 1 this function is strongly increasing 
with Ar. For 7 — 5/3, for example, dMp/dt oc — Ar^, 
whereas for molecular hydrogen, 7 « 1.4, and dMp/dt oc 
— Ar*. For definitiveness, we parametrise the mass flow 
rate through the LI point as 



dMp 
dt 



Ar 

Tp 



Mp 

tdyn 



if Ar > 



(7) 



where tdyn = \fr^jGMp is the (internal) dynamical time 
of the planet. Note that when the planet fills its Roche 
lobe, tdyn ~ 1/Slfe(a) = ^jGM,Jafi , and we can re- write 
equation [7] as 
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dMp _i _3/2 Mp 



Ar 



(8) 



where ao.i = a/(0.1 AU). Since the headhne factor in 
equation [8] is so large, one obtains very high Roche lobe 
overflow rates even for modest values of Ar. For ex- 
ample, even at Ar = O.Olrp the overflow rate is about 
two order of magnitude larger than the mean accre- 
tion rate onto embedd ed protostars, M ~ 10~* Mq yr"'^ 
(|Hartmann et alJlfOQ^ '). 

We found through experimenting with different 
power-law indexes in this dependence that the results are 
very insensitive to the exact form of the mass loss rate 
dependence on Ar. This result is well known from stellar 
binary evolutionary calculations, and will be re-derived 
for the problem at hand below. Physically, since the mass 
transfer rate is a strong function of Ar, it turns out that 
the process is stable only if |Ar| may be maintained very 
small, which is equivalent to requiring Vp ~ rn- In this 
situation the mass loss rate self-adjusts to satis fy the 
angular momentum loss or gain constraints (cf iRitteij 
1 19881 ). Different mass loss formulations only yield slightly 
different values of Ar at which the required equilibrium 
value of dMp/dt is achieved. 

On the other hand, in cases where the mass trans- 
fer leads to rn shrinking faster than rp, or r^ increas- 
ing faster than rn, a runaway mass transfer occurs. One 
well known example of this in stellar binaries is the 
case where the secondary is more massive than the pri- 
mary, in which case the separation of the components 
shrinks rather than in creases as the mass is transferred 
l|Pringle fc Wadelll985^ . In this case the disruption of the 
companion is nearly dynamical, and therefore the exact 
form of the transfer rate is again not very important. 



2.3 Roche lobe overflow: LI or both LI and L2? 

We follow here the usual assumption that the planet's 
sur face is determined by the equipotential surfaces (e.g., 
tof iFrank et al.ll2002D . We also assume that the planet 
and the star are on circular orbits around the centre 
of mass of the system. The dimensionless potential in 
a frame corotating with the centre of mass of the star- 
planet system has two saddle points, LI and L2. The 
points lie on the line connecting the star and the planet, 
with LI positioned between the planet and the star, and 
L2 further away "behind" the planet. 

In stellar binaries, the overflow of the secondary's 
Roche lobe always proceeds via the Lagrangian point LI 
since it is much closer to the secondary. However, for the 
planet-star system, the mass ratio, q = Mp/M,, is very 
small, and thus the overflow may proceed via both LI and 
L 2. We quantify th is following the treatment given in §4.1 
of lGu et al.1 (|2003l ). The difference between distances of 
the LI and L2 points {D\ and D2, respectively) from the 
centre of the planet are 

2rja 



AD = D2 



(9) 



where r? ^ rn /a = (Mp/SM.)^/^ < 1 (see eq. 63 in 
iGu et aL|[2003l ). 

We see that the point L2 is distance AD further away 



from the centre of the planet than LI. Now, if the atmo- 
sphere's scaleheight hp were infinitely small, the overflow 
would only proceed via the LI point (if the planet over- 
flow maintains the Vp = rn ~ Di condition). However, 
if the atmosphere is extended enough, that is hp ^ AD, 
then the outflow should proceed through both points LI 

and L2. 

Equation (66) of lCu et al] lj2003h shows that 



AD 



8 GMp 



(10) 



Our models will characteristically assume that Tcff ~ 10'^ 



K (see ^3.41 below), and the planet-star separation at 
which the Roche lobe overflow occurs is often a ~ 0.1 
AU. For these fiducial parameters. 



hp 
AD 



0.45 



Tcft 
10» K 



O.IAU 



Mp 



(11) 



This shows that for massive young protoplanets, Mp ~ 
lOMj, the overflow should indeed proceed mainly via the 
LI point, but as the planet's mass drops below IMj both 
LI and L2 can be letting the planet's atmosphere to es- 
cape. 

There is actually one further condition to check. In 
H2.2l we have shown that the overflow rate is a strongly in- 
creasing function of Ar. In the simulations below, we shall 
find relatively large Roche lobe overflow rates, requiring 
small but not infinitesimally small values for Ar/rp. The 
difference between locations LI and L2, AD, is also small 
in units of rp. For a large enough values of \dMp/dt\, it is 
possible that the required Roche lobe "overfilling", Ar, 
is larger than AD. We expect that if Ar ^ AD, the over- 
flow proceeds mainly via LI as argued above, but if Ar ^ 
AD, both LI and L2 will siphon gas away at comparable 
rates. To compare these, we note that for Mp = lOMj, 



AD 



TpTJ 



3/2 



0.06rp(Mp/10Afj)^/^ On the other 
hand, the planet mass loss of \dMp/dt\ = 10~* M© yr~^ 
requires Ar « 0.037rp. Thus we see that, at least for 
Mp = lOMj the overflow should proceed mainly through 
the LI point. 

In general, the ratio of the two length scales is 



Ar 
AD 



0.6 



\dMp/dt\ 
10-1 Mgyr- 



1/3 



1/2 



^ lOMj 



-5/6 



.(12) 



This shows that at the high mass loss rates the outflow 
proceeds via LI point as long as Mp ^ 5Mj; at lower 
planet masses the outflow may occur via both LI and L2 
points. 

We conclude that in the most typical parameter 
space sampled by our models below, the Roche lobe over- 
flow occurs via the LI point. However, for relatively low 
mass planets, Mp ^ IMj, the overflow should take place 
via both LI and L2. In this paper we make an approxi- 
mation that the outflow always occurs via LI. 



2.4 Definitions 

The time evolution of Ar is of primary interest to us 
here. Let us calculate the full time derivatives of the Hill 
(Roche) radius and the planet radius: 
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di 
and 

dlnvp 
dt ' 



dlnrn \ dlnA/p 
dlnM„ 



d In rp 
91nM„ 



dt 



dlnM„ 



+ 



a In 



dt 



+ 



dt 



d In rj 



(13) 



(14) 



respectively. Here, (d \nrp/d In Mp) is the adiabatic mass- 
radius exponent of the planet, discussed below, and 
{d\nrp/ dt) is the time derivative of planet's ra- 
dius due to thermal relaxation, e.g., radiative cooling 
of the planet. Similarly, (i91nr_fr/91nMp) is the rate of 
Hills radius change due to planet's mass changing, and 
{d In th /dt)j^ is the term due to external torques onto 
the planet in the absence of the mass loss. 

Since th = a{Mp/3Mt)^^^ , and in view of equation 

El 



dlnMp 



-2fa + 1/3 , 



(15) 



where we neglected the small term — {din M, /din Mp) = 
q = Mp/Mf <^ 1. Further, we introduce the following 
notations for brevity 

din r„ 



91nMp ' 

/ d In rp \ 
\ dt J 

din rn 
dt 



Mp=0 
1 



d In a 

~dr 



(16) 
(17) 
(18) 



The mass-radius relation for a polytropic planet is 
given by 

rp oc Mjf-" , (19) 

where n = 1/(7 — 1) is the polytrope's index, and 7 is 
the ratio of the specific heats for the planet's gas. It is 
noteworthy that the radius-mass relation for a polytrope 
strongly depends on 7, and that the polytropic cloud usu- 
ally (for large enough values of 7) expands as the mass is 
lost. For example, for 7 — 5/3, n = 3/2, and rp oc Mp 
for 7 = 7/5, n = 5/2, Vp oc M'^. 

Note that Tg can be related to the external torque. 
Text, on the planet. Here we anticipate that in our model 
situation the external torque is negative. Text < 0, push- 
ing the planet closer to the star. To relate Tc and Text, 
note that angular momentum conservation in the absence 
of mass loss yields 

Mp^ (GA/.a)^/" = -iTextl , (20) 

where (GAf.a)^''^ = Q aa^ is a specific angular momen- 
tum of the planet. From which we conclude that 



2lTc, 



MpO,aCl' 



< 



(21) 



Note also that Tc is negative for a planet con- 
tracting as the result of radiative cooling and it is 
positive for a planet expanding with time at a con- 
stant mass. The latter situation may occur when the 
planet's int erior is heated by the solid core accr etion 
luminosity (|Pollack et al.l Il996l : iNavakshinI l2011al '). or 



due t o dissipation o f tides induce d by the parent star 
(e.g.. IGu et al.ll2003l) . irradiation (|Burrows et al.l I2OO0I : 
iBaraffe et al.ll2003l ). 



2.5 Stability of mass transfer 

In accord with our plan to progress from simple to more 
complex we first study a non-cooling planet without ex- 
ternal torques, i.e., 1/tc = 1/tc — 0, which has Rh ~ Rp, 
i.e., transfers mass onto the star at some non-zero rate, 
— Mp^^ . As far as realistic planets in massive gas discs are 
concerned, this model situation is of an academic inter- 
est, as one needs external torques or internal evolution of 
the planet (e.g., planet swelling due to internal or tidal 
heating) to fill its Roche lobe, but it will be seen to pro- 
vide important hints to the results of the more realistic 
numerical experiments later on. 

Following the discussion in H2.2I the dominant factor 
in the mass loss rate dependence is the difference rp — rn- 
Write 

dlnM„ 



dt 



'(p{Ar) 



(22) 



where > is a positive and monotonically increasing 
function of its argument (note that dMp/dt < 0). We 
have explicitly neglected here the much slower depen- 
dence of the mass transfer rate on other variables (such as 
the photosphere's temperature and the density in equa- 
tion [S]). Within this approximation, we have. 



d^ In Mp 
dt2 



\ dt dt I 



(23) 



where cj>' = d(j>/dAr. Since at the onset of the mass trans- 
fer Tp'^TH, 

d^ In M„ .1 / dlnr-p d In th 

-<p Tp 



dt^ '"K dt dt J ' 

Referring now to eQuations ll3|[T^ and definitions 1151 \Wl 
we find that 



d^ In Mp _ . , dlnAfp 
= -0 rp (Cp - C-ff) 



(25) 



dt^ T-rV^K 

It is now apparent that the stability of the mass 
transfer rate is determined by the sign of (^p — C^h since 
d(j}' > 0. If this term is negative, the mass loss rate in- 
creases exponentially with time. Thus, the mass loss is 
unstable if 



Cp - Ch = Cp - ^ + 2/a < 



(26) 



where we used equation 1151 This relation demonstrates 
that stability of the mass loss process is decided by the 
mass-radius relation for the planet, and by the angular 
momentum transfer between the planet and the star (the 
factor fa)- Physically, this is because the planet moves 
outward as the result of the mass transfer, which in- 
creases rn- However, the planet may also increase in size 
due to mass loss. When the condition (|26p is satisfied, rp 
increases faster than rn as the planet siphons away its 
mass. In this case, once the planet filled its Roche lobe 
(Hill's radius), the mass loss process is a runaway one. 

The opposite situation, (p — (h > 0, is stable. In the 
presently considered case of zero external torques and no 
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cooling, the mass loss rate drops exponentially with time. 
The physics of this is simple - the planet moves outward 
and eventually it becomes smaller than its Roche lobe. 
In a more realistic situation where there are non-zero 
external torques, a steady-state situation can be set up, 
where the planet moves radially at just the right rate to 
maintain the balance ~ rn. 

Clearly, the stability of the mass transfer is strongly 
dependent both on the details of the evolution of the 
planet's angular momentum (described by the parame- 
ter fa) and by the evolution of the planetary radius in 
response to mass loss (described by the parameter ^p). 



2.6 Quasi equilibrium mass loss 

Let us now turn our attention to more realistic situations 
where the external torque and the planet's radius internal 
evolution are not negligible. In this case eQuation l24l leads 
to 



1 

H 



1 



(27) 



-lnM, = -0.,[(C.-CH)^ 

Since this is a linear differential equation, our conclusions 
about the stability of mass transfer do not change. If 
(p — < 0, the mass transfer is unstable as the mass 
loss rate runs away exponentially. 

However, if (p — (h > 0, there is a quasi steady state 
solution of the above equation with d^lnMp/dt^ « 0, 
where the mass loss rate is locked at a unique value given 

by 



(--- 



(28) 



This result, exc ept for differe nt notations, is the same as 
equation (16) of lRitteJ (| 19881 ). a well known result in the 
binary mass transfer studies. 

We note that if the planet is cooling and there is 
no internal energy source, then Tc < 0, e.g., the planet 
contracts with time (eguation llTp . The external torque is 
also negative by assumption, therefore for a contracting 
planet, 



Me, 



Mp 



Cp-Ch 



(29) 



equilibrium mass loss solution is possible only if |tc| > 
|re|. In the opposite case the planet contrac ts too quickly 
to be tidally disrupted (|Navakshinll201l"bl ): despite mi- 
grating in closer to the star, the mass transfer rate shuts 
down with time because Vp decreases faster than rn- 

On the contrary, if there are evolutionary reasons for 
the planet to swell with time, such as stellar irradiation 
or a powerful internal energy release by the solid core, 
then the thermal relaxation term in eg nation 128 1 amplifies 
rather than dumps the mass loss rate. In fact, mass loss 
from the planet may occur in this case even if there are 
no external torques on the planet. 

We can now justify our claim (also well known from 
the stellar binary mass transfer studies) in the end of ^2.2\ 
that the equilibrium mass transfer rate does not depend 
on the exact form of the mass transfer rate function, 0. As 
long as that dependence is monotonic and a strong one. 



the difference Sr = Vp — th adjusts to "the right value" 
which yields the equilibrium mass transfer rate given by 
equation 1281 The mass transfer rate in this sense is not 
a "driving" variable of the problem; it is one determined 
by the balance of the torques on the planet. 

Curiously, in the case of the external torques strongly 
dominating over the thermal relaxation of the planet, 
I Tel <C \tc\, the mass transfer rate is further insensitive 
even to the planet's mass, Mp. Consulting equation 1211 
we find that equation 1281 reads in this case. 



2|Te, 



That is, in this case, the mass loss rate depends only on 
the external torque and location of the planet, a, but not 
the mass of the planet. 

2.7 Radial migration due to mass loss 

Having understood the conditions under which the mass 
transfer rate can be steady, we now turn to the conse- 
quences of this for the radial migration of the planet. 
Adding the external torque on the planet to equation |3] 
we have 



d In a d In Mp 



1 

+ — . 



(31) 



dt ■' dt 

When the mass loss rate is in the steady state regime 
(equation I28|l . this becomes 



din a 
dt 



Cp - 1/3 



+ 



2fa 



To Cp - 1/3 + 2fa Tc Cp - 1/3 + 2/< 



(32) 



We remind ourselves that Tc is defined to be negative 
(equation I21[). that is, the the external torques are de- 
fined to drive the planet in rather than out. In contrast, 
Tc could be positive or negative. The remarkable property 
of equation[32]is that, depending on the different terms in 
the equation, the planet may migrate either in or out de- 
spite being pushed inward by the external torque, e.g., the 
outer disc. Furthermore, the relative magnitude of these 
terms may vary as the planet loses mass or evolves inter- 
nally. We shall later see that due to this, the same planet 
may change the direction of migration several times. 

2.7.1 Negligible cooling: inward or outward migration? 

When the thermal evolution of the planet is insignificant 
during the tidal disruption process, so that jTc| 3> Itc], 
we can neglect the cooling term in equation 1321 

din a 1 Cp — 1/3 



dt 



(33) 



Tc Cp - 1/3 + 2fa ■ 

Recall that in the steady-state mass transfer rate regime, 
the denominator of equation[33]is positive (cf. i)2.5p . Since 
Te < 0, we see that when the polytropic approximation 
for the planet is adequate, 

sign(^) = -sign (Cp- 1/3) . (34) 

That is, the planet migrates outward if Cp < 1/3, and 
inward otherwise. In particular, polytropic planets that 
expand as they lose mass (Cp < 0) always migrate out- 
ward. Physically, Cp controls the amount of mass that the 
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planet has to shed in response to being forced inward by 
the external torque. This response is meager for planets 
that shrink due to mass loss. In order to maintain the 
Tp « th equilibrium condition, such planets may simply 
lose mass as th shrinks. Thus external torques are able to 
push these planets in. In the opposite case, when < 0, 
applying the external torque to push the planet in causes 
it to lose so much mass that the outward torque of this 
material exceeds the originally applied torque, and the 
planet thus moves outward rather than inward. 



2.7.2 Rapidly cooling planets 

When thermal relaxation of the planet is not negligi- 
ble, e.g., \tc\ comparable to or much smaller than the 
external torque time scale, |rej, the last term in equa- 
tion [32] cannot be neglected. A casual look at the equa- 
tion [32] would suggest that a rapidly contracting planet, 
Tc < 0, \tc\ <^ |re|, migrates as prescribed by the last 
term of that equation. However, as remarked after equa- 
tion [29l the quasi-equilibrium mass loss is only possible 
if I Tel > I Tel or else the rapid thermal contraction of the 
planet shuts down the mass transfer entirely because rp 
becomes smaller than vh. 

Nevertheless, even in the quasi-equilibrium situation 
the second term in equation [32] may indeed exceed the 
first one and hence dictate the direction of the planet's 
migration. To give an example, consider our standard 
case of a planet with (^p = —1/3 (gas polytropic index 
7 = 5/3). Eg nation 132 1 reads for this case 



to the planet mass loss: 



din a 
dt 



1 



1 



Tc if a 



+ 



3fa 



1 



(35) 



Tc 3fa 

Therefore we see that for \tc\ somewhat larger than |re|, 
the planet may lose the mass in the quasi-steady regime 
and yet migrate outward if fa is sufficiently close to unity. 



2.7.3 Thermally expanding planets 

Let us now turn to the case of planets expanding as the 
result of thermal evolution, e.g., an inner energy source 
such as the solid core formation due to dust sedimenta- 
tion (|Navakshinl2011al ') . By definition, Tc > in this case, 
so the steady-state mass loss is possible for any value of 
To (provided i^p — i^j^ > 0) . The last term in equation [31] 
acts to push the planet outward; however, depending on 
the first term in that equation the migration direction 
can be inward or outward. 



_ 3__a_ 

dt ^ ROR 



^5(i?-i?dop) 



(36) 



where A = A/{flR)^, and A is the specific tidal torque, 
where 



A = - 















2 





R> a 



R<a. 



(37) 



In equation (|37p . f2 — ^ GMt / R^ is the angular veloc- 
ity at radius R, a is the radial position of the planet, 
q — Mp/Mt and p — R — a. This simplified form 
of the sp ecific torque is common l y used in literature 
(see, e.g.. iLin fc Papaloizou Il986l: Armitage fc Bonnelll 



I2OO2I : iLodato fc Clarkell2004l (Alexander et al.ll2006l ). We 

smooth the torque term for _R « a, where it would 
have a singularity (see equat ion I37II. We use the same 
smoothing prescription as in ISveTfc Clark3 (|l995l ) and 
[Lin fc Papaloizoul (|l986l ). i.e. for |7? — a| < max[ff, rj^], 
where H is the disc thickness and th = a{Mp/?,MtY^^ 
is the size of the Hill sphere (Roche lobe) of the planet. 

We shall note here that we use a constant a- viscosity 
parameterisation for the angular momentum transfer in 
the disc. Gravitoturbulent discs could in principle gener- 
ate an additional viscosity through gravitational t orque s 
(jLin fc Pringlel 1 19871 : iGammid |200]J : [Rice et al.l |2005| ). 
However, here we are interested in the inner few AU discs 
that are unlikely to be self-gravitating since the radiative 
cooling is very ine fficient in units of local dynamical tim e 
in such discs (e.g., lRafiko"^l2005l : [ciarke fc Lodatdl2009l ). 
Such small discs are also expected to be much less mas- 
sive than the central star and henc e the global model of 
mass transfer is als o very unlikely (jLodato fc RicdboOSl : 
ICossins et"al]|2009l ). 

The last term in equation [36] describes the mass 
deposition in the inner disc by the planet (note that 
Mp < 0). The mass deposition takes place into a nar- 
row ring centred onto the deposition radius i?dop, which 
is the circularisation radius of the planet's material lost 
through the LI point into the inner disc. The specific an- 
gular momentum of the gas at LI is Qa{a — th)^ , where 
Q,a = a/ GMt/a'-^, and thus the deposition radius is given 
by 



3 NUMERICAL METHODS 

3.1 Mass and angular momentum transfer in 
the disc 

We use the code of iLodato et al.l (|2009l ) for the evolution 
of an accretion disc in the presence of an embedded satel- 
lite, rescaled to the case of a protostar of mass A/* and 
a planet of mass Mp. The disc is described by a diffusive 
evolution model that includes the tidal torque term aris- 
ing from the planet and the mass deposition term, related 



i?dcp = a (1 - ^) . (38) 

Operationally, we find the radial zone inside of which 
i?dcp fall and deposit all of the mass lost by the planet 
there. Note that J?dep is a function of time. Depositing gas 
in a single zone does not appear to lead to any numerical 
problems as the gas spreads in both directions from that 
locations fairly quickly, resulting in a well behaved surface 
density profile. 
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3.2 Thermal disc equations 

The rest of disc equation s follow the standard vert ically- 
averaged approach (e.g.. IShakura fc Sunvaevlll973h with 
minimum necessary complexity. In particular, disc vis- 
cosity in the a-prescription for a gas-pressure-dominated 
disc is 



(39) 



where Cg = kT^/ fimH is the gas sound speed at the disc 
midplane, Tc is the midplane temperature, p is the di- 
mensionless mean molecular weight and mn is the mass 
of the proton. The disc vertical scale height H is found 
by solving for the vertical pressure balance, taking into 
account the radiation pressure (although for most cases 
considered here radiation pressure is generally negligible): 



GA'L 



R 



-pc 



(40) 



where pc and Pgas ~ PckTc/ fimn are the gas midplane 
density and pressure, respectively, and arad is the radia- 
tion energy density constant. The central gas density and 
the disc surface density are related by E = 2pcH . 

The midplane equilibrium temperature is solved as- 
suming the vertical energy balance for the disc: 



(41) 



9,(r.,)Ef7^ . ^ , 

where r = «:(pc, Tcq)pcH is the disc opacity for the opac- 
ity coefficient Hi{pc,Tcc^. For opacities varying slowly with 
temperature, this approach is entirely satisfactory, as 
the thermal equilibrium is established on the time scale 
fth = l/ofi, which is much shor ter than the disc viscous 
time (e.g., cf. lFrank et"alll2002l '). 

We modify this standard equation to account ap- 
proximately for the irradiation heating of the disc by the 
star, whose luminosity is set to L, = LQ(il'/*/M0)^. 
The incidence angle of the irradiation for simplicity is set 
to give cosi = 0.1 everywhere, although in principle this 
should be a function of position and the disc vertical scale 
height, H{R). Thus we write 



^irr - 

and 



/ L. cosi y/"* 



(42) 



(43) 



For protostellar discs, however, due to a strong tem- 
perature dependence of opacity coefficient k on tempera- 
ture (primarily) , the disc may be thermally unstable (e.g. 
iBeU fc Lin|[l99i ). In terms of equation 1431 it means that 
there may be more than one solution for a given E, and 
that some of the solutions may be thermally unstable. 
The time evolution of the disc in this case depends on 
its previous thermal state, and one needs to solve a time- 
dependent version of the energy equation rather than as- 
sume thermal equilibrium. 

While thermal instabilities are clearly relevant to the 
problem at hand, we would like to concentrate here on 
the "new" physics - the planet tidal disruption and the 
resulting inner disc refilling - as much as possible. We 
therefore pick a very simple and transparent numerical 



method to deal with thermal instabilities in the present 
paper. In this minima list approach, similar to that of 
iLodato fc Clarke! (|2004l ). we use a time-dependent energy 
balance equation with a radial advection of heat term: 



dT 
'dt 



T — T 



dT 
OR 



(44) 



In this equation, Toq is found from equation |l3] using the 
previous timestep value for the gas temperature, Tt, in 
the disc opacity calculation, k(pc, Tt). The radial velocity 
vr is given by the following: 



d 



Vr 



T,Ry^ OR 



[uJlR 



1/2] 



+ 2nRX 



(45) 



Our simple energy equation allows us to capture the ther- 
mal disc instability if it does occur without an excessive 
associated numerical cost. 



3.3 Radial motion of the planet 

In our approach, there are two ways in which the an- 
gular momentum is passed between the planet and the 
disc. First of these is due to gravitational torques (the 
second term on the right hand side of equation I36|l and 
is independent of the instantaneous mass loss rate by the 
planet. The term represents the back reaction of the disc 
on the planet's motio n and follows from a ngular momen- 
tum conservation (see lLodato et al]|2009l '): 



R^ \T.dR (46) 



•.at /J 



■JRin 

where the integral is taken over the whole disc surface and 
Qa is the angular velocity of the planet. The subscript 
"J" indicates that the term arises due to the planet-disc 
angular momentum exchange. 

The second channel for the angular momentum ex- 
change is only present when the planet loses mass to the 
disc. As explained in i)3.1l the disc gains material with the 
specific angular momentum of the LI point of the planet. 
The material is deposited at the respective circularisation 
radius, -Rdop (equation [38} . The planet's specific angular 
momentum thus increases at the rate 



= --Mp [n„a^ - Qa{a - rnf] .(47) 



Adding the two terms for the radial migration of the 
planet, we obtain 



din a 
dt 



VGALa 



Mp 2rH 



2 - 



XT.dR 



(48) 



Note that the planet mass loss rate appears directly only 
in the second term on the right hand side, but it is also 
indirectly present in the first one. As the planet fills in 
the inner disc with gas, it changes the disc surface density 
profile, and therefore the first term on the right as well. 
This system of equations is thus quite non-linear. 
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3.4 Planet contraction 

In this paper we are only interested in tiie planets dis- 
rupted inside the inner 1 AU of the s tar. Su ch disruptions 
were termed "hot" by iNavakshinI (|2011bl ) due to their 
possible physical links to the "hot" jupiters and lower- 
mass close-in exoplanets. Only giant planets in which 
molecular hydrogen has been dissociated into atomic or 
ioni sed H are dense enough to venture into this region 
(see iNavakshinlboid ) . 

Our isolated constant m a ss plan et model is based on 
the toy model of INavakshinI (|2011bl ) who assumed that 
the young plane t follows a Hayashi-li ke track during its 
early evolution (jCraboske et al.lll975l ). The model gives 
the radius of the planet, rp{t), at time t, as 



4{t) = 



1 + Arlt 



(49) 



where r2 is the planet's radius at the time of the second 
collapse, and A = 2'i7raBT^g / (GMp) , where the effective 
temperature of the planet, Teft is a free parameter that is 
fixed for a given simul ation. To a c commo date a range of 
possible cooling rates, INavakshinI (|2011bl l explored mod- 
els of contracting planets with three values of the effective 
temperature, 500, 1000 and 2000 K, but we test only the 
mid range model here. 

The initial radius of the planet, r2, right after molec- 
ular hydrogen dissociation is also uncertain, with re- 
sults depending on dust opacity, internal energy liber- 
ation by a possibly present solid core, planet rotation, 
and the detailed e quation of state at high densities (cf. 
iMarlev et allbOOTl ). Following iNavakshin (2011bj), we es- 
timate the initial (post-collapse) radius of the planet 
by the energy conservation argument. We assume that 
the change in the gravitational potential energy of the 
planet, ~ —GMp/2rp, as the planet radius collapses, is 
mainly used to dissociate H2 molecules and ionise hydro- 
gen atoms. Molecular hydrogen is completely dissociated 
after the collapse, whereas the fraction of fully ionised 
hydrogen, ^ ^ 1, is a free parameter of the model 
that encapsulates the uncertainties pointed out above . 
Going through this simple logic (cf. INavakshinI [2011bl ). 
we obtain for the initial radius of the planet. 



r2 



GMpTUH 

D + 2Xi£ 



Mj 1 + 2X,£/D 



(50) 



where D = 4.5 eV and £ — 13.6 are the dissociation 
energy of hydrogen molecules and the ionisation energy of 
hydrogen atom, respectively. The lowest starting density 
models correspond to Xi = 0, whereas Xi = 1 give the 
densest possible models. 

In equal ion [49] the time t is counted from the time of 
the second collapse, which we take to be the start time of 
the simulations presented below. This simple model does 
not take into account the electron degeneracy pressure. 
For low density start models, electrons b ecomes partially 
degen erate only after t ^ 10® yrs (e.g., iGraboske et al.l 
119751 ). We do not run our models for that long in the 
present paper. For the high density start models, we over- 
estimate the density of the planet (since the electron de- 
generacy pressure is neglected), but this becomes impor- 
tant only in the innermost regions of the disc, R ^ 0.02 



AU, where the planets are likely to be swallowed by the 
star anyway. We also neglect the importance of the possi- 
ble solid core inside the planet. We plan to address these 
issues in future publications. 

In practice, we start our calculations with the planet 
of a given mass Mp right after the second hydrodynamical 
collapse. The initial planet's radius, rp{0) = r2, is given 
by equation [SD] When the mass loss from the planet com- 
mences (cf. !j2]2}, the planet's radius is advanced accord- 
ing to equation 1141 The right hand side of that equation 
is the sum of the term due to the planet's mass loss and 
the radiative cooling term calculated at a constant planet 
mass. For our toy cooling model, we have, 



dlnvp 
dt 



3 ^"^"^ ' 



(51) 



where A = 2'iTTcyBT^s/{GMp{tf ), i.e., calculated with 
the evolved Mp{t). For reference, 



1.8x10" yr 



fJ^f fO-Ol ril^V.(52) 



V IQMj 



I- 



3.5 Parameter /„ and the gap in the disc 

Our analytical study showed that the outcome of the tidal 
disruption of a planet strongly depends on the parame- 
ter fa that describes the fraction of the specific angu- 
lar momentum lost through LI but then regained by the 
planet due to gravitational torques from the inner disc 
(cf. equation |3}. In numerical experiments below, these 
torques are calculated self-consistently, so fa is neither 
needed nor introduced. But it may be quite beneficial to 
know what an effective value of fa is for a given numerical 
situation to be able to appeal to the analytical insights. 

The simplest situation here is when the planet carves 
out a deep gap in the disc that stops any significant gas 
flows across the gap. The outer disc can then be viewed 
as a source of an external torque onto the planet, and 
the situation is thus quite analogous to that in a semi- 
detached stellar binary system. In a binary system, in the 
conservative mass and angular momentum exchange sce- 
nario, the disc transfers mass in and angular momentum 
out, giving the excess speciflc angular momentum to the 
secondary. In our problem, if the star accretes gas with 
the specific angular momentum of the last circular orbit 
around the star (which we set for simplicity to R = R,), 
V GMfRf , and the planet recovers the rest of the specific 
angular momentum lost through the LI point, then 

- {R,/af/' . 



fa 



(53) 



If a > 7?., then fa ^ 1. 

Now, when the gap is partially or completely closed, 
gas is able to flow past the planet, taking with it the 
excess angular momentum as well. While details clearly 
depend on how deep the gap is, we expect that /„ be- 
comes smaller as the gap closes, and in fact goes to zero 
in the limit of no ga£; 

According to iTakeuchi et al.l (|l996l l. the planet 
opens an annual gap in the disc if Mp exceeds 



1/2 



(54) 
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As H/R ~ 0.1 — 0.2, the estimated gap opening mass 
varies from a sub- Jupiter to a few Jupiter masses depend- 
ing on the a-parameter of the disc. Obviously, various de- 
grees of stellar irradiation and different opacity laws can 
affect the gap opening through increasing or decreasing 
the geometrical aspect ratio H/R. 



4 TESTS OF PLANETARY DESTRUCTION 

We now present a series of controlled numerical experi- 
ments. To focus our attention on the planet only, we do 
not yet introduce an outer disc. Instead, we express the 
effects of the outer disc by a fixed torque, 

f ^) = - , (55) 

\ at J cxt "^c 

where = —250 yrs for the tests below. Except for the 
extreme mass ratio, Mp/M, , the setting of the problem is 
exactly analogous to that of mass transfer in stellar bina- 
ries, where the external torque might be due to the grav- 
itational radiation or magnetic breaking. Such a setup 
facilitates a straight forward comparison with the analyt- 
ical theory for the planet migration and Roche lobe over- 
flow developed earlier. However, we shall find "new" be- 
haviour compared with stellar binaries motivated by the 
disc gap closing (something that never happens for stellar 
binaries). More self-consistent simulations, those without 
an external torque and where the outer disc torques the 
planet in, are to follow in iJS] 

In this section, the planet initial mass is Mp — lOAfj, 
the initial planet-star separation is ao = 0.25 AU, the 
initial planet's radius, Rp = 0.021AU « 447?j, obtained 
from equation [50] with Xi = 0. The effective temperature 
of the planet is set to TcS — 1000 K. We do not allow the 
planet's mass to fall below A/corc = 5M0, assuming that 
there is a solid core of that mass inside the planet. This 
particular assumption only affects the results at the very 
end of the tidal disruption of the planet. We use Thomp- 
son scattering opacity for the disc for simplicity (we shall 
use a more realistic one for the full disc simulations later 
on). The inner radius of the computational domain, as- 
sumed to equal the stellar radius, is Rin = 0.01 AU, and 
the outer radius is -Rout = 4.2 AU (which is sufficient for 
these tests as we are only interested in the inner regions). 
The number of grid zones is Nr — 500, unless otherwise 
stated. 

Table 1 shows the list of the simulations presented 
in this paper, most important parameter values, respec- 
tive figure numbers and also the most notable behaviour 
patterns of the disc or the planet. 

4.1 Simulation Extl 

The first run presented, labelled Extl, is for a disc vis- 
cosity parameter a — 0.02 and the planet radius-mass 
exponent of — —1/3. 

4.1.1 Evolution of the planet 

Figure [T] shows the most interesting characteristics of this 
simulation. The solid line in the top panel shows the 



planet's location, a, as a function of time, in units of 
AU. The dashed curve shows the planet mass in units 
of lOMj. The solid curve in the middle panel shows the 
radius of the planet, Rp, whereas the dashed curve shows 
the Hill's radius. The bottom panel of Fig. [T] presents 
three mass loss/gain rates: the accretion rate onto the 
star, the absolute value of the planet mass loss rate, \Mp\, 
and the analytical prediction for that mass loss rate given 
by equation 1281 shown with the solid, the red dotted and 
the blue dashed curves, respectively. 

The general outcome of the simulation is that the 
planet spirals in for about 200 years with zero mass loss, 
until it fills its Roche radius and starts losing mass. The 
mass lost by the planet is transferred into the disc inte- 
rior to it. A quasi-steady state is established in which the 
mass loss by the planet is exactly matched by the mass 
accreted onto the surface of the star. The inner disc trans- 
fers its excess angular momentum back to the planet, 
producing an outward directed torque on the latter. This 
reverses the direction of planet's migration. Losing mass, 
it is migrating outward, as predicted. The quasi-steady 
state mass loss cannot last forever as the planet's mass is 
finite, and eventually the mass loss runs away, destroying 
most of the gas component of the planet rapidly (see the 
spike at around 660 yrs). Only ~ 0.12Mj of the gaseous 
planet survive the disruption in this simple model (see 
igTS] below). 

The blue dashed curve in the lower panel of Fig. [T] 
shows the predicted mass loss rate given by [2H1 with fa 
given by equation [S3] We observe that the analytical pre- 
diction provides an excellent description of the numeri- 
cal results starting from the time when the planet fills 
its Roche lobe, e.g., from time t « 200 yrs to t « 550 
yrs. During this time the mass loss rate from the planet 
is matched by the accretion rate of the star. The sec- 
ond marker of this quasi-equilibrium is the fact that the 
rate of planet's expansion/contraction due to mass loss 
and radiative cooling is exactly matched by the corre- 
sponding change in the Hill's radius, so that the two are 
approximately equal (cf. the middle panel of Figure [T} . 
Note that the radius of the planet increases from t ~ 200 
to « 400 yrs, and then decreases until the sharp spike. 
This change in behaviour is explained by the fact that 
initially the planet's radius evolution is controlled by the 
mass loss, with cooling being negligible. Thus the planet 
initially expands. After ~ 400 yrs, when the planet loses 
about half of its original mass, the radiative cooling be- 
comes dominant in controlling the planet's radius evolu- 
tion. The planet then contracts until the accretion rate 
spike. 

The quasi-equilibrium state of the planet-star system 
is analogous to that of a semi-detached stellar binary sys- 
tem where the lower mass secondary transfers its mass to 
the primary. However, after t ~ 550 yrs the planet's mass 
loss rate stops decreasing with time, flattens off and then 
increases. The accretion rate onto the star follows the 
trend, but lagging behind the planet's mass loss. As we 
shall see below, this is not a time-delay effect but rather 
is a signature of a very important effect that does not oc- 
cur in stellar binaries: the closing of the gap in the disc, 
due to which the inner disc material may flow past the 
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Figure 1. Evolution of a massive gas planet tidally disrupted 
near its parent star. The top panel shows the radial loca- 
tion and planet's mass versus time, as labelled. The middle 
panel shows the planet's radius, rp, and its Hills radius, ru- 
The bottom panel shows the accretion rate onto the star, the 
planet mass loss rate and the theoretical prediction for that 
rate (equation I28I I with the solid black, red dotted and blue 
dashed curves, respectively. The theoretical prediction works 
exceedingly well in the region where the planet fills its Roche 
lobe and the gap in the disc is opened (see i|4.1.1ll . 

planet. As the result of this, the planet suffers a nearly 
catastrophic disruption at f « 660 yrs. 

4-1.2 Evolution of the disc 

Figure [2] shows several snapshots of the disc surface den- 
sity versus radius covering the most interesting moments 
of this numerical experiment. The curves are marked by 
the respective time t shown right next to them. As stated 
earlier, initially the surface density profile is everywhere 
zero. At t « 200 yrs (see the middle panel of Figure [l]), 
the planet eventually fills its Roche lobe and starts trans- 
ferring mass into the inner disc, filling it out. The first 
curve shown in Figure [2] corresponds to time t = 198 
yrs (black solid curve), when the planet disruption just 
begins. The planet at that moment is at a = 0.12 AU. 
The material lost through the LI point is circularised and 
is deposited at 7?dep ~ 0.06 AU (note the break in the 
black curve in the Figure), and then spreads viscously 
in both directions. At small radii the material eventually 
starts accreting onto the star, whereas at larger radii the 
torques from the planet prevent it from spreading out- 
wards further. During this time the mass loss rate from 
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Figure 2. The disc surface density for several different times, 
as marked on the Figure. The evolution of the disc can be 
summarised as following: (a) the inner disc is inflated with 
material from the disrupted planet; (b) the disc pushes the 
planet outward; (c) the gap is partially closed, leading to a 
runaway in planet's mass loss rate; (d) the disc overflows the 
remaining low mass planet. The evolution after that is that of 
a "standard" accretion flow unperturbed by the planet. 

the planet exceeds the accretion rate onto the star (the 
rising part of the red dotted curve at f ~ 200 yrs in the 
bottom panel of Figure [ij . The inner disc is thus being 
filled in by the material lost from the planet. 

The swelling of the inner disc cannot go on forever, 
however. When the mass of the disc becomes significant 
enough [t = 211 and t = 219 yrs curves in Figure [2]), 
the gravitational torque of the inner disc starts to affect 
the radial migration of the planet. The inner disc torque 
in fact exceeds that of the externally applied torque and 
the planet reverses its direction of migration. Migrating 
outward, it is now able to settle into the self-regulated 
quasi-steady regime discussed earlier. 

This amiable quasi-steady state evolution of the 
planet and the inner disc continues until time t « 500 — 
600 yrs or so. By that time, the mass of the planet drops 
to only ~ 3Mj. This severely undercuts the ability of 
the planet to keep the gap open (cf. equation [54]). This is 
why the gap starts to be partially filled by gas, and the 
material from the inner disc starts to leak into the outer 
one (cf. the t — 358 yrs curve in Figure [2]). By t — 643 
yrs curve, the gap is seriously compromised. The last of 
the curves, t = 971 yrs, shows that in the end the gap is 
completely erased and the disc spreads viscously outward 
even more. By that time the planet has lost most of its 
mass and the disc is hardly aware of its presence (even 
though it consists entirely from the material formerly be- 
longing to the planet in this simulation). 

4-1.3 The late spike/final destruction of the planet 

The partial gap closing at t ~ 600 yrs has catastrophic 
consequences for the planet. As the inner disc material 
diffuses outward through the gap, past the planet's orbit, 
less of the inner disc torque is passed on to the planet. In 
terms of the analytical parameterisation of ij2l the frac- 
tion of the angular momentum passed back to the planet, 
fa, drops. Since (p — C,H ~ 2(/a — 1/3) for 7 = 5/3 planet. 
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the mass transfer rate from the planet becomes unstable 
when fa drops below 1/3. The mass loss rate runs away 
exponentially at time t ~ 660 yrs. Physically, the planet's 
radius increases more rapidly due to the mass loss than 
the Hill's radius could increase since the outward radial 
migration is compromised by the leaking gap. As can be 
seen from the bottom panel of figure (TJ during this expo- 
nential runaway, the protostellar accretion rate does not 
manage to keep up with the rate at which the matter 
is deposited in it. This is not surprising as some of the 
matter flows outward through the gap. 

The runaway planet destruction via Roche lobe over- 
flow is bound to end in one of the two ways. One is a 
complete unbinding of the gaseous envelope that would 
leave only a rocky core (which is present for this simu- 
lation by assumption). On the other hand, if the cooling 
time of the planet becomes very short at low masses, 
the planet may go into a "runaway contraction" phase 
instead, where it contracts faster than the Hill's radius 
expands. This second outcome is indeed what happens in 
simulation Extl. The end mass of the planet is 0.12Mj, 
including the assumed 5 of the rocky core. 

We emphasise that the final mass of the planet 
strongly depends on its internal structure, and so is very 
model dependent. In our simple constant Tcb — 1000 K 
cooling prescription, contraction time scale becomes very 
short at low planet's masses. This shuts off the mass loss 
at the end of simulation Extl. A more realistic cooling 
model would probably lead to longer cooling times. On 
the other hand, the solid core may have created its own 
dense gas "atmosphere" strongly bound to the core, in 
analogy to the nucleated c ore instability model for gi- 



Simulation Ext2 



ant p lanet formati on (e.g., Mizun ol ll980l : IPoUack et al.l 



tP 

1 19961 ). as argued bv lNavakshinI (2011a). This atmosphere 
is not modelled here, and may remain bound to the core. 
Clearly, better models of the internal structure of the 
planets must be used to address the issue of the post- 
disruption mass. 



4.2 Importance of gap opening: simulation 
Ext2 

As we have seen in simulation Extl, when the gap closes, 
fa drops since the matter and angular momentum flow 
from the inner disc into the outer disc, past the planet. 
This reduces the rate of the outward migration of the 
planet and may destabilise the mass loss by the planet, 
leading to its runaway disruption. In view of the gap 
opening criterion (equation 154 1) . then, the outcome of the 
planetary tidal disruption depends on the disc viscosity, 
the planet mass, and the outer disc torque. The latter 
controls the accretion rate in the inner disc (through the 
equilibrium mass loss condition in equation I30p which in 
its turn regulates the geometrical disc aspect ratio, H/R. 

To explore this issue, we run one further simulation, 
labelled Ext2, entirely identical to Extl except that the 
viscosity parameter a is increased to a — 0.1. Figure [3] 
shows the resulting time evolution of the planet in the 
same format as Figure [1] The outcome of the higher a 
simulation is markedly different from that of its lower a 
counterpart. There is no quasi-steady state plateau in the 
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Figure 3. Same as Fig. [T] but for simulation Ext2 with a 
higher disc viscosity parameter, a = 0.1. Due to this, the disc 
is hotter, and the planet is unable to keep the gap opened. 
As the result, the material from the inner disc is able to flow 
past the planet. The reduced angular momentum feedback on 
the planet destabilises the mass transfer process. Most of the 
planet's gaseous inventory is completely disrupted on a very 
short time scale. 



accretion rate corresponding to an equilibrium tidal evap- 
oration of the planet. Analysis of the disc surface density 
evolution (not shown here for brevity) for simulation Ext2 
shows that there is only a depression in the surface den- 
sity profile rather than a deep gap. As explained earlier, 
a partially opened gap lowers the effective fa- The planet 
does not move outwards quickly enough, and the expan- 
sion of the planet (recall that Rp oc Mp^^^) increases the 
mass loss rate further, until the planet is destroyed in a 
runaway fashion. 



4.3 Importance of the mass— radius relation: 
simulation Ext3 

As a final numerical experiment in a series with the fixed 
external torque timescale, we present simulation ExtS, in 
which the mass-radius relation for the planet is steeper, 
rp oc Mp^ . That is, ("p = —1 rather than —1/3 for simu- 
lations Extl and Ext2. Figure |4] shows the outcome of 
this simulation, demonstrating that there are actually 
two mass loss episodes and both are unstable. 

The stability of the mass transfer depends (see ^2.Q\ 
on the difference C,p — C,h = "if a — 4/3 in this case. If the 
gap is impenetrable to the gas from the inner disc, and 
the material lost by the planet is instantaneously gained 
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Figure 4. Same as Fig. [T] but for simulation Ext3. Here the 
mass-radius relation of the planet is steeper, Vp oc Mp^. The 
mass-loss from the planet fuels too fast an expansion of rp, 
so that the difference rp — ru "runs away", and so does the 
planet mass loss rate. The planet loses most of its mass on 
a very short time scale in the first very luminous burst. The 
remaining gas envelope is unbound in the second disruption 
wrhen the planet migrates inward even closer. 



by the star, then /a = 1 — \J R*/a, as explained in i|3.5l 
The planet is first disrupted at R ~ 0.11 AU in this 
experiment, so /„ ^ 0.7, and Cp ^ Ch ~ 2fa — 4/3 ^ 0.1, 
i.e., small but marginally positive. However, there is a 
small lag between the mass loss rate of the planet and 
the accretion rate, which reduces the effective value of 
fa below the above (theoretically maximum) estimate. 
This is probably the reason for the mass transfer being 
unstable in this simulation. 

The first disruption episode is not final in this simu- 
lation. As the planet loses all but ~ O.SA/j of its gaseous 
mass (which corresponds to only 3% of its initial mass), 
the radiative cooling accelerates, and rp shrinks rapidly. 
At the same time, the planet migrates outward, so the 
Hill's radius increases. This allows the remaining Saturn- 
mass planet to shut down the mass loss for a while. How- 
ever, the inner disc drains onto the star, while the "exter- 
nal torque" in our fixed Te model continues to push the 
planet inward. So the planet resumes its inward migra- 
tion. By the time t = 800 yrs it finds itself ai R = 0.06 
AU, where it fills its Roche lobe for the second time. A 
second outburst, physically but not quantitatively sim- 
ilar to the first one, results. The remaining reservoir of 
the gas is lost. The "naked" solid core then continues to 
migrate inward as prescribed by the fixed Te model. 



5 DISRUPTION OF PLANETS IN 
REALISTIC DISCS 

5.1 Parameters and initial conditions 

We now consider a more realistic set of calculations in 
which we do not introduce an artificial external torque. 
Instead, we start with an accretion disc filling in the com- 
putational domain. The torque parametrized in SjUby the 
timescale Tc is now self-consistently computed from the 
time-dependent disc structure. The disc is initialised with 
a surface dens ity profile commonly used in the literat ure 
(e.g.. Matsu vama et al.ll2003l : [Alexander et al] 120061 ') in 
studies of protostellar disc evolution: 



So(ii) = ^ 






■ R ■ 


exp 


. Ro. 



(56) 



where Ro is the disc length-scale, set at i?o = 5 AU un- 
less stated otherwise, and Am is a normalisation con- 
stant chosen so that the disc contains a given initial mass 
Md = 27r RdRT,o{R). In hindsight, the exact shape 

of the initial surface density profile does not appear very 
important, as long as one samples all the interesting pa- 
rameter space by varying other parameters of the prob- 
lem, such as the disc mass Md and the v i scosity param- 
eter a. We use the opacity of IZhu et al.l (|2009l ) instead 
of Thompson electron scattering opacity, used in tests 
Extl-Ext3. The computational domain is modelled by 
Nr = 500 to 700 radial zones. We also use = 0.02 AU 
for the inner boundary of the disc, instead of 0.01 AU as 
used in tests Extl-Ext3. The outer boundary of the com- 
putational domain is set to Rout = 25 AU for the tests 
below, or otherwise explicitly stated. We use a reflecting 
boundary condition at i?out. This choice is of no partic- 
ular importance for the results as there is little gas mass 
at large R and the most interesting effects take place at 
R ^ 1 AU. In addition, the viscous time at i?out is much 
longer than the duration of the simulations below. The 
initial mass of the planet is set at Mp — lOMj. 



5.2 Simulation Siml 

Simulation Siml is done with the following parameters: 
viscosity parameter a = 0.002, Md = 20Mj, Xi = 0, 
Toff = 1000 K, Cp = -1/3. The initial location of the 
planet, ao = 0.13 AU, is chosen to be such that the Roche 
radius rn is just slightly larger than the planet radius, 
rp. This is done for computational expediency. Tests Sim3 
and Sim4 below are started with larger value for ao. 

Figure [5] shows the resulting evolution of the planet 
and the accretion rate onto the star. In brief, the time 
evolution of the system can be described as the following 
stages: (i) From t — to t ~ 100 yrs. The evacuation of 
the disc interior of the planet; (ii) t ~ 100 to 400 yrs. 
An inward migration of the planet and then commence- 
ment of the Roche lobe overflow, which re-fills the inner 
disc; (iii) t ~ 400 to 2600 yrs; the quasi steady state out- 
ward migration of the planet as the result of the inner 
disc torque exceeding that of the outer disc torque. The 
planet loses most of its mass during this time, evolving 
from lOMj to only about IMj. (iv) t 2660 yrs. The gap 
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Figure 5. Same as Fig.[T]but for simulation Siml tiiat models 
the disc torques self-consistently. 

is partially closed. This reduces the inner disc torque and 
destabilises the quasi-equilibrium condition [r-p ~ rn)- 
The planet mass loss runs away since expands faster 
than m could increase due to the outward migration. 
This leads to the "final destruction" spike in the planet 
mass loss rate (the dotted curve in the bottom panel of 
figure [5| . The planet loses almost all of its gaseous mass 
in this particular simulation, (v) t > 2660 yrs, the "post- 
planet" disc evolution. With the planet now all but de- 
stroyed (the remaining solid core is dynamically unim- 
portant for the disc), the disc evolves independently of 
the planet. 

We now discuss this evolution in greater depth. Fig- 
ure[6]shows the disc surface density profile at key epochs. 
Panel (a) concentrates on the earlier times, when the 
planet is still massive enough to keep the gap opened. 
The solid curve shows the initial E(_R). The next curve 
at t = 81 yrs shows a greatly reduced surface density 
in the disc interior to the location of the planet. It also 
demonstrates that the planet pushes the outer disc back, 
creating a bump in S(i?) just behind the gap. After the 
inner disc disappears (accretes onto the star), the outer 
disc pushes the planet closer to the star. Roche lobe over- 
flow sets in, and the inner disc is refilled by the planet's 
material. The t = 712 curve in Figure [6^ shows that it 
is filled almost to the same values of E as at t = 0. The 
planet now migrates outward, and the gap becomes less 
wide since the planet's mass is reduced. This description 
(and Figure [S^) covers simulation stages (i) to (iii), as 
defined above. 

Figure [6]d shows what happens when the planet's 



mass is too small to keep the gap opened. The gap is 
significantly compromised by the time of the black solid 
curve, t = 2654 yrs. The gap is completely closed, and the 
planet's mass is dispersed and assimilated into the inner 
disc in the next ~ 20 yrs (see the t — 2665 yrs curve). 
The mass released by the planet streams both inwards - 
onto the star - and outwards (note the spike in the red 
t = 2665 yrs curve in the figure). 

Figure [7| zooms in onto the time interval around the 
"final destruction" of the planet, e.g., around the gap 
closing, showing the planet's radius and Hills radius in 
the top panel, and the accretion rate onto the star (solid) 
and the planet's mass loss rate (dotted) in the bottom 
panel. The planet's destruction in this simulation ends 
when only ~ 0.1 of gas remains bound to the 5M(g 
solid core we assumed in this simulation. Given our toy 
cooling model (fixed Tcff), this amount of mass cools very 
quickly, allowing the radius to contract very fast and 
avoid a complete unbinding of all of the gaseous mass. It 
is clear that this outcome is completely dependent on the 
cooling model, which becomes unrealistic at this point. 
Nevertheless, this demonstrates an important point. If 
radiative cooling becomes sufficiently rapid as the planet 
loses mass, the planet (the gas atmosphere of the solid 
core, actually) may contract below the Hill's radius be- 
fore the evaporation of all of the gas. The end result of 
this is then a solid core plus a gas envelope planet. 

We see from Figure [7| that most of the material lost 
by the planet in the "last disruption" outburst is accreted 
by the star within half a century or so (see also the curve 
t = 2831 yrs in fig. [Sb)- The disc then evolves indepen- 
dently of the planet, accreting onto the star. The sharp 
transitions in E(7?) in Figure [S)d at late times {t = 6425 
yrs, for example) are caused by the strong opacity jumps 
where Hydrogen becomes partially ionised. 



5.3 Siml.X: More compact start models 

So far we explored only the low density start models, in- 
troduced by setting the parameter Xi — (cf. equation 
I50p . This parameter stands for the fraction of the ionised 
hydrogen atoms in the planet right after the second hy- 
drodynamical collapse. In our toy model for the planet 
this parameter controls the radius of the planet right af- 
ter H2 dissociation. The larger the initial radius of the 
planet, the farther away from the star can the planet be 
disrupted. We now explore how the results depend on the 
value of Xi , keeping the other parameters of the simula- 
tion fixed at same values as in Siml (cf. Table 1). We 
label the simulations Siml.X, where X = 10 x Xi. 

Figure [5] shows the outcome of simulation Siml. 4, 
with Xi = 0.4. Since the initial radius of the planet is 
less than one third that of simulation Siml (compare 
the solid curves in the middle panels of Figs. [S]and[S|, 
the planet is able to migrate farther in before it fills its 
Roche lobe. The planet is located at 7i ~ 0.035 AU at 
time t « 2300 yrs when it fills its Roche lobe. In con- 
trast to simulation Siml, the mass transfer rate is not 
quasi-steady in this case. The whole of the planet's en- 
velope is tidally disrupted in a short Roche lobe overflow 
event. The main destabilising factor is the decreased frac- 
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Figure 6. The disc surface density profile for several snapshots from the simulation Siml. The times of the snapshots are labelled 
on the figure next to the respective curves, (a) The early evolution of the system. The t = initial condition neglects the presence 
of the planet, but the planet opens a gap in the disc very quickly. The inner disc then empties out onto the star. The outer disc 
then "pushes" the planet closer in where it fills its Roche lobe and starts transferring mass to the inner disc. The inner disc then 
refills with the planet material. Note that the gap becomes narrower as the planet loses mass, (b) Later evolution of the disc, 
centred around the spike at t = 2660 yrs (see Figs. [5] and [7]l. The planet is rapidly dispersed when the mass loss rate runs away 
due to the partially closed gap. An ionisation front propagates outward, but stalls quickly. After the giant planet is destroyed, the 
disc accretes onto the star as if there were no planet. 



tional amount of the angular momentum that the planet 
receives back when its mass is accrete d by the proto- 
star. In particular, /a = 1 — y'RhJa ~ 0.24 for the 
present simulation, whereas for Siml fa ~ 0.55 due to the 
more distant disruption location. For (^p = —1/3, we have 
i^p — (^H ~ 2{fa — 1/3). Therefore, for simulation Siml, 
— Ch > 0, whereas for Sim 1.4 — C^h < 0, explaining 
the vastly different outcomes of these simulations. 

The mass deposited in the inner disc drives the 
planet outwards to about 0.2 AU. During this short 
planet Roche lobe overflow spike our model for the mass 
deposition is probably somewhat unrealistic, as the mass 
may also overflow through the L2 point, in the outer disc, 
at these very high outflow rates. 

Figure |9] shows the stellar accretion rate evolution 
near the tidal disruption point for simulations Siml, 
Siml. 2, Siml. 3 and Siml. 4. The rise times of the light- 
curves are between a year to ten, but the more compact 
the planet is (the larger Xi), the steeper the rising and 
the falling parts of the curves. This is natural, as the 
more compact models are disrupted closer in. The rise 
time of the light-curves is about the inner disc viscous 
time, which becomes shorter at smaller values for deposi- 
tion radius, R, scaling as cx R^^^ (at a fixed a parameter 
and H/R). 



5.4 Simulation Sim2: higher disc mass 

We now present simulation Sim2, which is identical to 
Siml except for a higher initial accretion disc mass, Md = 
50Mj. Figure [TOl shows the main results of this numerical 
experiment. Since the disc accretion rate is higher, the 
torques acting on the planet are higher too. Therefore 
the time scales for planet migration in this simulation are 
shorter than in Siml. In terms of our analytical approach 
to the problem, the external disc torque timescale, |re|, 
is shorter in Sim2 than in Siml. 



As before, the inner disc is first emptied onto the 
star. The planet fills its Roche lobe at around t = 100 
yrs, starting a very large accretion event. Note that the 
maximum accretion rate reached in this simulation at the 
first Roche lobe overfiow event is nearly 10~* M© yr~^, 
almost an order of magnitude higher than in simulation 
Siml, which is explainable by the larger disc torques onto 
the planet. 

During the quasi-steady Roche lobe overflow stage, 
the planet migrates outward, as before. Interestingly, by 
the time it reaches IMj mass, the planet has migrated 
outward to i? ~ 0.2 AU, farther than it did in simu- 
lation Siml. This can be understood from equation 1321 
The two terms on the right hand side of that equation 
have different signs, with the last term (due to radiative 
cooling) causing inward migration. This term remains 
roughly the same in Siml and Sim2 as Vp is almost the 
same in the two simulations. However, the external disc 
torque is much larger in Sim2. Since the net result of the 
outer disc torque is outward migration in this case, the 
planet migrates outward faster in Sim2 than it does in 
Siml. 

We would like to pause here for longer to clarify 
the physical interpretation of the faster radial migration 
in Sim2 compared with Siml better. At first this result 
seems to contradict the common sense: "A more mas- 
sive outer disc should provide a larger spin down torque, 
so the planet should migrate outward slower or even mi- 
grate inward in Sim2 as compared to Siml" . However, the 
equilibrium mass loss conditions imply that the planet 
must lose mass more rapidly in Sim2 precisely because 
the outer disc torque is larger. This means that the in- 
ner disc in Sim2 is more massive than it is in Siml, and 
it produces a larger outward torque on the planet. To 
be in the equilibrium situation, the outward and inward 
disc torques on the planet must work in such a way as 
to push it steadily outward at a rate proportional to the 



© 2011 RAS, MNRAS 000.[Tlf23l 



16 Sergei Nayakshin and Giuseppe Lodato 




10'°^ : 

2640 2660 2680 2700 2720 
time [yrs] 



Figure 7. A zoom in on the evolution of the planet's radius, 
and the Hills radius (top panel), and the planet mass loss rate 
and the star's accretion rate (the bottom panel) for simulation 
Siml. 
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Figure 8. Same as Fig. \5\ but for simulation Siml. 4 which 
differs from Siml by a higher hydrogen ionisation parameter, 
Xi = 0.4. The planet is thus about 3 times more compact, thus 
being disrupted closer in to the protostar. The disruption is 
abrupt and non steady because C^p — (^h < for this simulation 
(see i|5.3l for detail). 



outward disc torque (cf. equation [33]) . So physically, a 
larger outer disc torque causes a stronger response from 
the planet through the inner disc feeding and that is why 
the planet migrates outward even faster. 

As the result of this more effective outward migra- 
tion, the planet does not actually suffers a nearly catas- 
trophic disruption when it reaches M « lAfj. At the 
start of the second planet mass loss spike (time t « 400 
yrs in fig. llOf) . the planet moves outward quickly. Instead 
of suffering a runaway mass loss as in Siml, a shutdown of 
the mass loss occurs at this moment in Sim2. The planet 
continues to be pushed outward by the inflated inner disc, 
reaching i? ~ 0.3 AU at t ~ 500 yrs. By that time the in- 
ner disc drains sufficiently onto the star, and the inward 
migration of the planet resumes. Despite some contrac- 
tion of the planet's radius with time, the Hills radius 
becomes smaller than again at time t ~ 800 yrs. This 
time Roche lobe overflow is terminal for the gaseous part 
of the planet, and it is completely destroyed in the third 
Roche lobe overflow event very quickly. The evolution of 
the disc after that is similar to Siml. 

5.5 Thermal instability in the gap (Sim3) 

To explore the importance of initial conditions, we set 
up simulation Sim3, which has same initial disc mass as 
Sim2, Md = 50Mj, but higher disc viscosity parameter. 



a = 0.01, and a different disc radial cutoff Ro = 7 AU (see 
equation I56p . The outer boundary of the disc is enlarged 
to i?out = 40 AU. We also place the planet at larger 
starting separation, ao = 1 AU rather than 0.13 AU in 
Siml, Sim2. 

The results of this test are shown in Figure 1111 If 
we were to average out the shorter period (tens of years) 
variations in the curves, the gross evolution of the sys- 
tem is not that different from that found previously. Once 
again, the inner disc drains onto the star whilst the outer 
disc is "dammed up" behind the planet. The accretion 
rate onto the star then plummets to very small values 
until the planet is pushed to R ~ 0.11 AU where it 
overfills its Roche lobe. The material siphoned from the 
planet feeds the inner disc, restarting the accretion onto 
the star. Neglecting the oscillatory pattern for now, the 
planet loses mass and moves first outward and then in- 
ward until Mp m 1.5Mj, all the time satisfying rp ~ rn- 
When the gap is partially filled, there is a very large 
outburst that unbinds the residual gaseous mass of the 
planet. After this, the gap in the disc shuts completely, 
and the planet is of no importance for the disc. 

However, on shorter time scales, there is far more 
variability in Sim3 than in any presented before. Figure 
[12] zooms in onto two selected time intervals from simula- 
tion Sim3, showing just the mass loss and accretion rates. 
An oscillatory behaviour is seen already ~ 70 years af- 
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Figure 9. Comparison of the protostellar accretion rates, M», 
around the first Roche lobe overflow event for simulations Siml 
and Siml.X. The simulations differ from one another only by 
the initial ionised hydrogen fraction, Xi = 0, 0.2, 0.3 and 0.4, 
for the dashed, solid, dotted and dash-dot curves, respectively. 
Note that the more compact the planet is (the larger Xi is) , the 
shorter the rise time of the outburst and the larger the peak 
value of Mt . The curves were arbitrarily translated along the 
time coordinate to appear on the same figure. 
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Figure 10. Same as Fig.[5]but for simulation Sim2. Note that 
the timescales are shorter and the mass loss rates arc larger 
in Sim2 than they arc in Siml. 
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Figure 11. Same as Fig. llOl but for simulation Sim3, in which 
the planet is inserted at a = 1 AU instead of a = 0.13 AU. The 
oscillations in the mass loss/gain rates in the bottom panel are 
duo to thermal ionisation instability in the gap, as explained 



ter the the first Roche lobe overflow (cf. the top panel of 
fig. ll2|l . Limit-cycle variations persist throughout most of 
the time preceding the "final destruction spike" (cf. the 
bottom panel of fig. I12p. 

Figure [T51 shows the disc central (midplane) temper- 
ature, surface density, and the effective temperature, TcS, 
at two times, t = 1171 yrs, and t = 1195 yrs. The first of 
these times (solid black curves in Figure I13|l is close to 
the third peak in the stellar accretion rate curve in the 
top panel of fig. 1121 whereas the second set of curves (red 
dashed curves in fig. 1131 respectively) corresponds to the 
next trough in the accretion rate. 

Concentrating on the top panel of fig. 1131 we note 
that the two curves differ by a factor of ~ two in the 
inner disc, and are almost identical behind the gap re- 
gion. However, the gap region is significantly different at 
these two selected times. In particular, the disc central 
temperature Tc ~ 2 x 10** K at the peak of the accretion 
rate cycle, and only Tc ~ 800 K at the minimum of the 
cycle. 

These profound changes in the disc temperature have 
similarly profound consequences for planet evolution and 
the inner disc, as that is strongly linked to the planet's 
migration and mass loss rate. The middle panel of fig. [TS] 
shows that while the gap is only partially opened at the 
"high-T" state of the gap region, it gets fully opened at 
the "low-T" state, when viscosity of the material in the 
gap drops in response to the temperature variation. 
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Figure 12. Same as the bottom panel of Fig. llll but zooming 
onto two interesting time periods in the disc/planet evolution. 
In the top panel, the beginning of the outburst is shown, when 
the Roehe lobe is overfilled for the first time. The bottom 
panel of the figure corresponds to the period around the "final" 
destruction of the gaseous envelope of the planet. 



However, what causes the temperature drop in the 
gap region? To address that, note that the planet mass 
loss rate actually dips much faster than the stellar ac- 
cretion rate after any one of the peaks visible in the top 
panel of ll2l This implies that the planet has been pushed 
outward "too far" by the inner disc torques, and it con- 
tinues to be pushed even further out for some ~ 10 yrs 
after the peak. The inner disc eventually runs out of ma- 
terial, as it is accreted onto the star. The outer disc then 
over-powers the torques from the inner one and the planet 
migrates inward again, restarting the mass loss and re- 
filling the inner disc for another cycle. Thus, while the 
gap changes are pronounced, they appear to be driven by 
the behaviour of the inner disc, which makes the planet 
to switch between the full-gap to a partial-opened gap 
states. 



6 PLANET-DISC THERMAL 
INSTABILITIES 

During the outburst state, the inner disc becomes 
very hot, with most of hydrogen atoms ionised. One 
may therefore expect that the a-parameter may be 
high er, in accord with estimates from dwarf nova sys- 
tems iKing et al.l (|2007|). and Magneto- Ro tational Insta- 
bility simulations ( Balbus &: Hawlevll 19981 ). Furthermore, 
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Figure 13. Disc behaviour through one of the instability cy- 
cles for simulation Sim3. Top; disc central temperature; Mid- 
dle: disc surface density; Bottom: disc effective temperature. 
The curves are for two selected times, t = 1171 yrs (solid black 
curves), and t = 1195 yrs (dashed red curves). These corre- 
spond to a peak and a trough in the star's accretion rate (cf. 



modeling of the observations of FU Ori outbursts by 
IZhu et al.l (|2007l ) suggested that a ~ 0.02-0.2. We there- 
fore present now simulation Sim4, which is identical to 
Sim3 except a = 0.04. 

Figure [14] shows the outcome of this simulation, 
which contains both similarities and also very significant 
differences from all of the simulations presented earlier. 
Similarly to the previously studied cases, the inner disc 
first empties onto the star, while the outer disc pushes the 
planet in. However, unlike any of the previous cases, the 
first accretion outburst episode onto the star is not driven 
by the planet filling its Roche lobe and transferring its 
mass inward. Instead, at time t ~ 179 yrs, thermal disc 
instability is triggered. At that time, the planet is situ- 
ated at i? « 0.2 AU and does not fill its Roche lobe. 

To analyse the disc behaviour in greater detail, Fig- 
ure[T5]shows the disc structure just before the instability 
is triggered (black solid curve), during the initial luminos- 
ity rise (red dashed curve) and just after the rise (blue 
dotted). One notices that the instability is actually trig- 
gered at the outer banked-up edge of the gap, behind 
the planet, at i? ~ 0.3 AU. Before the instability, the 



© 2011 RAS, MNRAS 000.[Tlf23l 



Fu Ori outbursts and planet-star mass transfer 19 



/rds/tag/data1/sn85/FU_ORI/Mp10_Md50_alpha0.04_XiO_a1/ 





200 400 



Figure 14. Same as Fig. 1141 but for simulation Sim4, which 
is identical to Sim3 except that Q-parameter is increased to 
a = 0.04. The resulting behaviour change is significant: the 
thermal disc instability is triggered before the planet over- 
flows its Roche lobe. The resulting disc/planet behaviour is 
somewhat similar to that found by Lodato & Clarke (2004). 



disc is relatively cold everywhere, but especially so in the 
gap. As the instability flares up, the disc behind the gap 
heats up and the increased pressure of the gas is able to 
partially close the gap. The gas thus rushes inward, ac- 
cumulating in the inner disc, sending that region of the 
disc on the hotter, upper stabl e branch of the "S-curve" 
(|BeU fc Linlll994l : lLodato fc Cla rke 200i), and leading to 
a massive rapid-rise outburst. 

The outburst in t he simulation Sim4 is in a qual- 
itative agreement with iLodato fc Clarkj (|2004l ). except 
here it is supplemented by the mass loss from the planet. 
Note that in this high a (and thus high accretion rate) 
simulation, the mass loss from the planet is a relatively 
minor detail. Indeed, when the mass loss from the planet 
sets in (t « 300 yrs), it exceeds the prior value of the ac- 
cretion rate onto the star only in the flrst ~ 50 yrs or so. 
Together with Sim3, this simulation clearly demonstrates 
that non-linear coupling between the hydrogen ionisation 
instability in the disc and the planet-disc mass exchange 
may be important in determining a variability pattern in 
the accretion rate onto the star. 



7 EMISSION REGION SIZE 

Figure [TS] shows the disc central temperature (top panel) 
and the disc effective temperature (bottom panel) for sev- 
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Figure 15. The disc structure at times t = 176.4 yrs, 179.2 
and 179.5 yrs for simulation Sim4, shown with the black solid, 
red dashed and blue dotted curves, respectively. Note that 
thermal hydrogen ionisation instability is triggered just behind 
the gap, and the instabilit y then propagates in, as in models 
of llLodato fc Clarkell2004l) . 



eral representative times for the simulation Sim4. Along 
with those curves, we also s how the best fit power- l aws o f 
the inner hot disc model of lEisner fc Hillenbrandl (|201ll ) 
to their NIR Keck interferometer data for the three well 
known FU Ori objects. These observations have spatial 
resolution as good as ~ 0.05 AU, and thus serve as direct 
and sensitive probes of the spectral disc model fits. 

The effective temperature profiles of our discs in the 
inner regi on are not that dissimilar f rom the best fit 
models of lEisner fc HiUenbrandl (|201lh . except for the 
presence of the gap. This qualitative agreement is en- 
couraging, given that we did not try to fine tune our 
models to satisfy the spectral constraints. This suggests 
that our models may be relevant to the observed FU 
Ori sources, although it is desirable to attempt more de- 
tailed spectral fits to the observations in the spirit of 
lEisner fc Hillenbrandl j^OlJ). This is however beyond the 
scope of the present paper. 
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Figure 16. The disc central temperature (Top panel) and the 
effective temperature (Bottom panel) versus radius for simula- 
tion Sim4, at three selected times, t = 1036, 1248 and 1584 yrs, 
shown with the thick solid black, red dashed and blue dotted 
curves, respectively. The three black solid power-laws repre- 
sent the best fitting effective temperature profiles inferred by 
Eisner and Hillenbrand (2011) for the three FU Ori sources. 

8 DISCUSSION 

8.1 Main results of this paper 

In this paper we presented calculations of the tidal dis- 
ruption of massive gaseous planets in massive accretion 
discs of young protostars. In the analytical part of the pa- 
per we derived conditions under which the mass and 
angular momentum exchange between the planet and the 
disc is in a quasi steady state (cf. text around equation 
I26p . and the corresponding equilibrium mass loss rate 
from the planet. Qualitatively, these conditions require 
(a) the planet to be massive enough to keep the gap in 
the disc opened; and (b) the mass-radius relation for the 
planet to be such that the planet contracts as it loses 
mass or at least expands only as a weak power-law of the 
planet's mass, Mp. 

We then formulated a numerical ID approach to 
modeling the planet-disc mass and angular momentum 
exchange (?SJ. Utilising a toy cooling model for a young 
planet, and a fixed external torque driving the planet in. 



we found that our numerical model reproduces the ana- 
lytically expected results well O, including the quasi- 
steady state planet's disruption (Roche lobe overflow), 
and the eventual closing of the gap and a runaway planet 
disruption when the planet's mass becomes too low. We 
then proceeded to several tests (labelled Siml to Sim4 in 
Table 1) in which the planet-disc torques are calculated 
self-consistently. 

The most robust result of our calculations is that if 
there is a massive young planet in the inner disc, then 
the accretion rate onto the protostar must be strongly 
variable on time scales ranging from a few years to a few 
hundred or even thousands of years. There are simply 
too many non linear physical processes occurring in the 
planet-disc system for the accretion rate to remain static 
over long time scales. 

First of all, as well known from previous literature 
l|Rice et al.ll2003l ) a massive planet opens a deep gap in 
the disc, blocking the flow of matter from the exterior 
disc to the inner one. Such a blockade of the protostar 
by the young planet may be the reason why the aver- 
age accretion rates of p rotostars appear to b e lower than 
expected theoretically (|Dunham et al.ir2010l ). 

We found two ways in which the blockade is ended, 
and both lead to a rapid rise accretion outburst that may 
be relevant for the observations of FU Ori outbursts. The 
first of these, encountered in the relatively high a simula- 
tio n Sim4 (a = 0.04, cf. jj5] and figure [T^ . is the same as 
the lLodato fc Clarke! (|2004 ) model for FU Ori outbursts. 
In this case thermal (hydrogen ionisation) instability sets 
in at the outer edge of the gap in the material piled up 
behind the planet. Heating up, the outer disc is able to 
close the gap, overfiow the planet and fuel an outside-out 
accretion outburst. 

The second way to kick-start accretion onto the star 
was found in all of our lower a simulations (Siml to Sim3 
in Table 1). In these cases the planet is pushed inward 
sufficiently fast to fill its Roche lobe before the hydro- 
gen ionisation instability at the outer gap takes place. 
The planet then loses mass through its LI point. The 
material circularises inward of the planet, and spreads 
viscously. As gas accretes onto the protostar, the accre- 
tion rate rises by orders of magnitude. What happens 
next depends on the internal properties of the planet. If 
conditions for the quasi-steady state mass transfer are 
satisfied, then the planet migrates radially (usually in 
the outward direction), so that the Hills radius is exactly 
equal to the planet's radius. The resulting outburst light 
curve may then resemb le a step funct ion, similarly to the 
FU Ori outburst itself (|Clarke et aLllioOS) . 

During this quasi-steady state planet-star mass 
transfer, the role of the planet completely reverses. While 
before the outburst the planet prevented the material 
from the outer disc from reaching the star, now it is the 
donor for the star. Also, like in a stellar binary system, 
the planet takes the excess angular momentum of the in- 
ner disc away and does not allow the material to spread 
outward viscously. This latter effect keeps the accretion 
rate onto the star at a higher value for longer than it 
would have been if the gas could flow past the planet. 

If the steady-state conditions for mass transfer are 
violated, then an even more dramatic outburst may oc- 
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cur. The planet is then tidally disrupted on time scale of 
a few to few tens of years. The material disrupted off the 
planet then spreads viscously in both directions, and the 
accretion rate onto the star drops rapidly as there is no 
planet to bank up the inner disc in the outward direction. 

On the top of this generic behaviour, we found vari- 
ability on smaller time scales caused by a non-linear con- 
nection between the planet migration and mass loss and 
hydrogen ionisation instability. For example, in simula- 
tion Sim3 ( ^5.5p , only the gap region of the disc is unsta- 
ble and switches between the ionised and the non-ionised 
states, modulating the accretion rate onto the protostar. 

These results do depend quite sensitively on the in- 
ternal structure of the planet, as we described in i\3A\ 
In particular, if planets are smaller after H2 dissociation 
than we assumed in most of our simulations here, then 
they may not be disrupted at all. In this case the planets 
are driven all the way into the star or the magnetospheric 
cavi ty. This correspon ds to the "high density start" mod- 
els (|Navakshinll2011bl ) parametrised by Xi = 1. There- 
fore, better modeling of the internal relaxation of the 
growing planets, including possible massive solid cores 
and the luminosity released by these, is absolutely neces- 
sary to detail predictions of our model further. 



8.2 Where do the inner planets come from? 

It is possible that every young protostar goes through 
as many as 10-20 FU Or i events during its growth 
(|Hartmann fc Kenvon|[l99^ ) . Recent mod eling of the em- 
bedde d phase of low mass proto stars by iDunham et al.l 
(|2010l ) also suggests that large amplitude variability in 
the accretion rate is widespread, and that stars may ac- 
tually acquire most of its mass in bursts. If our model for 
FU Ori outbursts is to apply to the observations, many 
giant planets are required per star. Is this reasonable? 

The answer strongly depends on which planet forma- 
tion theory is accepted. In the standard Core A ccretion 
(CA) model (|Pollack et al.lll996l : lAlibert et al] [2005;) of 
planet formation, giant planets are formed "late", e.g., 
after a few Myrs. These planets are also expected to 
be rather compact (|Marlev et al.ll2007l '). Tidal disruption 
of giant CA planets is therefore unlikely, and having as 
many as ten of such per star would also seem implausible. 

The gravita tional disc instability model for planet 
formation (e.g., iBosd [l997[ ) was thought to form plan- 
ets only at distances greater than 50-100 AU (e.g., 
iRafikovij gOOS'') due to inefficient c ooling at closer dis- 
tances jGammic 200ll : iMaver et a"l. 2004; Rafikov 2005); 
iMeru fc Bat6l201ll ). It may therefore appear that plan- 



ets born at these large distances would be of no rel- 
evance to this paper. However, working in the con- 
text of proto stellar acc r etion rather than planet forma- 
tion, IVorobvov fc Basul (|2005l . I2OO6I I showed that their 
clumps migrated inward rapidly, reaching all the way to 
their inner computational domain boundary of 10 AU. 
These results ha ve been supported by their more re- 
cent simulations (|Vorobvov fc Basul I2OIOI ) a nd by simu- 
lations of a number of independent authors (iBolev et al.l 
2OIOI: iBolev fc DurisenI l20ld: ICha fc Navakshinl I2OIII 



Machida et al.ll201ll : iBaruteau et al.ll201ll : lMichael et all 



I2OIII ) . IBaruteau et al.l (|201ll ) provide an in depth analy- 
sis of the rapid migration of gas clumps in self-gravitating 
discs and conclude that these clumps migrate via type I 
regime. These authors find that massive planets migrate 
inward faster than a disc gap could be opened. A recently 
submitted manuscript by Zhu et al (2011b; private com- 
munication) also shows many examples of massive clumps 
migrating inward rapidly although the authors find that 
some of these clumps become substantially more massive 
than a giant planet. 

Note that these results only apply to the outer 
gravito-turbulent discs; in the inner non self-gravitating 
disc regions studied here, planet migration must slow 
down considerably just because planets we consider are 
1-3 orders of magnitude more massive than the total mass 
of these inner discs. Standard gap opening criteria should 
app ly then. 

iBolev et all (|2010l ); iNavakshinI |20lO) argued that 
tidal disruption of gaseous envelopes of the giant proto- 
planets may be a plausible way to form all kinds of plan- 
ets in what was called the "Tidal Dow nsizing Hypothesis" 
for planet formation. More recently, INavakshinI (|2011bl ) 
have shown that proto-planets in which molecular hy- 
drogen is dissociated are dense enough to undergo tidal 
disruption inside the inner ~ 0.1 AU from the star. 

If these ideas are to be relevant to FU Ori outbursts 
and episodic accretion of stars, one must accept that 
most early formed planets are completely destroyed and 
swallowed by their parent stars. This is required both 
from the fact that many ou tbursts are needed per star 
(|Hartmann fc KenvonjligQel ). 



9 CONCLUSIONS 

In this paper we studied the planet and the disc evolu- 
tion in cases when the young massive gas giant planet 
overfills its Roche lobe and transfer its mass back into 
the disc. Due to significant uncertainties in the internal 
structure of young planets, large parameter space for the 
problem, and still missing physics, our study cannot yet 
give definitive answers on what exactly happens with the 
planets. Instead, this work is to be viewed as one estab- 
lishing an analytical and numerical framework for further 
studies of the planet-disc-star mass exchange. 

It appears that Roche lobe overfiow of young giant 
gas planets holds a potential promise as an explanation 
for the enigmatic FU Ori outbursts of protostars. We 
note in passing that a subset of our models (such as 
Sim3) may be relevant to shorter period and lo wer ac- 
cretion rate variable protostars such as EXors jHerbigj 
Il989l ; ISiciha-Aguilar et al.|[200i : iLorenzetti et al.[|2009l ). 
If FU Ori outbursts and episodic accretion onto proto- 
stars are indeed connected to tidal disruptions of young 
planets then the implication would be that many young 
giant planets are tidally destroyed or perhaps swallowed 
whole by their parent proto-stars. 

Future work should include better treatment of the 
planet's internal structure, especially a possible presence 
of truly massive solid cores, add type I migration (ne- 
glected here) , the possibility for the overflow through the 
L2 point (neglected here), and aim to model the "end" 
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Name " 


(Mj) 


a" 




(AU) 




Figs. B 


Comments'' 


Extl 


- 


0.02 





0.25 


-1/3 


[Tl2l 


Reference case 


Ext2 


- 


0.1 





0.25 


-1/3 


El 


Gap closed 


Ext2 


- 


0.02 





0.25 


-1 


H 


Unstable mass transfer 


Siml 


20 


0.002 





0.13 


-1/3 


[5i6im 


Reference case 


Siml.X 


20 


0.002 


X 


0.13 


-1/3 


[8l9l 


closer-in disruptions 


Sim2 


50 


0.002 





0.13 


-1/3 


M 


Temporary shutdown of mass transfer 


Sim3 


50 


0.01 





1 


-1/3 


[mini [13] 


Thermal instability in the gap 


Sim4 


50 


0.04 





1 


-1/3 


M 


Disc instability before mass loss 



Table 1. List of simulations performed and main parameters. Notes: 

" Simulation IDs. Those starting with "Ext" are simulations with an imposed external torque ([|4]l. Those with "Sim" are full 
disc simulations where the disc-planet torques are calculated self-consistently. 
^ Initial disc mass for the full disc simulations. 
Cf-parameter of the simulation. 

The fraction of ionised hydrogen used to determine the initial planet radius (eq. I50I I . 

The starting location of the planet, in AU. 
^ The index of the planet's mass-radius relation (equation [T6J . 
s Corresponding figures in the paper. 
^ Notable behaviour in the simulation. 



result when the gas disc is presumably dispersed by disc 
phot o- e vap or at ion . 
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